!========================================================================================================
 program timeaverage
! 
! Compute the time-averaged profiles of a few quantities
! The program needs the file time_profiles.dat that can be produced by running
! make_timeprofile.csh
!
! Run it as follows 
! ./timeaverage 'output_dir' time1 time2 nz
! where
! - 'output_dir' is the directory where 'time_profiles.dat' was stored (by default, make_timeprofiles.csh
!   creates this file in the directory where simulatio outputs are stored)
! - time1 and time 2 are the initial and final non-dimensional times over which the time-averages must be computed
! - nz is the number of vertical grid points
!========================================================================================================
	implicit none 
	integer :: i,j,k,nz
	real (kind=8) :: time,t1,t2
	real (kind=8), dimension(:), allocatable :: zgrid,vrms,T,eta,Phi,comp,ux,uz,RaH,alpha,kap, &
	                                            vrms_tav,T_tav,eta_tav,Phi_tav,comp_tav,ux_tav, &
						    uz_tav,RaH_tav,alpha_tav,kap_tav
        character (len=100) :: dirout, arg						    
!=========================================================================================================

        do i = 1, iargc()
           call getarg(i, arg)
	   if(i.eq.1) dirout=arg
	   if(i.eq.2) read(arg,*) t1
	   if(i.eq.3) read(arg,*) t2
	   if(i.eq.4) read(arg,*) nz
        end do

! Set the time interval over which the time-average must be computed
        nz=nz-2 
	allocate(zgrid(0:nz+1),vrms(0:nz+1),T(0:nz+1),eta(0:nz+1),Phi(0:nz+1),comp(0:nz+1),  &
	         ux(0:nz+1),uz(0:nz+1),RaH(0:nz+1),alpha(0:nz+1),kap(0:nz+1))
	allocate(vrms_tav(0:nz+1),T_tav(0:nz+1),eta_tav(0:nz+1),Phi_tav(0:nz+1),comp_tav(0:nz+1), &
	         ux_tav(0:nz+1),uz_tav(0:nz+1),RaH_tav(0:nz+1),alpha_tav(0:nz+1),kap_tav(0:nz+1))

	open(1,file=trim(dirout)//'time_profiles.dat')

        vrms_tav=0d0; T_tav=0d0; eta_tav=0d0; Phi_tav=0d0; comp_tav=0d0; ux_tav=0d0; uz_tav=0d0
   	RaH_tav=0d0; alpha_tav=0d0; kap_tav=0d0
	k=0
	do 
	   do j=0,nz+1
	     read(1,*,end=666) time,zgrid(j),vrms(j),T(j),eta(j),Phi(j),comp(j),ux(j),uz(j),RaH(j),alpha(j),kap(j)

	     if(time.ge.t1.and.time.le.t2) then
	       backspace(1)
	       if(j.eq.0) k=k+1
  	       read(1,*) time,zgrid(j),vrms(j),T(j),eta(j),Phi(j),comp(j),ux(j),uz(j),RaH(j),alpha(j),kap(j)
	       vrms_tav(j)=vrms_tav(j)+vrms(j)
	       T_tav(j)=T_tav(j)+T(j)
	       eta_tav(j)=eta_tav(j)+eta(j)
	       Phi_tav(j)=Phi_tav(j)+Phi(j)
	       comp_tav(j)=comp_tav(j)+comp(j)
	       ux_tav(j)=ux_tav(j)+ux(j)
	       uz_tav(j)=uz_tav(j)+uz(j)
	       RaH_tav(j)=RaH_tav(j)+RaH(j)
	       alpha_tav(j)=alpha_tav(j)+alpha(j)
	       kap_tav(j)=kap_tav(j)+kap(j)
	     endif

	   end do
	end do
	
666	vrms_tav=vrms_tav/k
        T_tav=T_tav/k
	eta_tav=eta_tav/k
	Phi_tav=Phi_tav/k
	comp_tav=comp_tav/k
	ux_tav=ux_tav/k
	uz_tav=uz_tav/k
	RaH_tav=RaH_tav/k
	alpha_tav=alpha_tav/k
	kap_tav=kap_tav/k

	open(2,file='timeave_profiles.dat')
	do j=0,nz+1
	  write(2,80) zgrid(j),vrms_tav(j),T_tav(j),eta_tav(j),Phi_tav(j),comp_tav(j), &
	              ux_tav(j),uz_tav(j),alpha_tav(j),kap_tav(j)
	end do
80 	format(20E20.8)


 end program timeaverage


 
